R Notebook (Earthquake)

Lo primero que hacemos es cargar todas las librerías que vamos a utilizar.

library(FSA)

## ## FSA v0.8.25. See citation('FSA') if used in publication.
## ## Run fishR() for related website and fishR('IFAR') for related book.

library(Rmisc)

## Loading required package: lattice

## Loading required package: plyr

library(psych)

## 
## Attaching package: 'psych'

## The following object is masked from 'package:FSA':
## 
##     headtail

library(boot)

## 
## Attaching package: 'boot'

## The following object is masked from 'package:psych':
## 
##     logit

## The following object is masked from 'package:lattice':
## 
##     melanoma

library(rcompanion)

## 
## Attaching package: 'rcompanion'

## The following object is masked from 'package:psych':
## 
##     phi

library(e1071)
library(fitdistrplus)

## Loading required package: MASS

## Loading required package: survival

## 
## Attaching package: 'survival'

## The following object is masked from 'package:boot':
## 
##     aml

## Loading required package: npsurv

## Loading required package: lsei

library(ggpubr)

## Loading required package: ggplot2

## 
## Attaching package: 'ggplot2'

## The following objects are masked from 'package:psych':
## 
##     %+%, alpha

## Loading required package: magrittr

## 
## Attaching package: 'ggpubr'

## The following object is masked from 'package:plyr':
## 
##     mutate

library("car")

## Loading required package: carData

## 
## Attaching package: 'car'

## The following object is masked from 'package:boot':
## 
##     logit

## The following object is masked from 'package:psych':
## 
##     logit

## The following object is masked from 'package:FSA':
## 
##     bootCase

library(nortest)
library(classInt)
library(gridExtra)
library(plotly)

## 
## Attaching package: 'plotly'

## The following object is masked from 'package:ggplot2':
## 
##     last_plot

## The following object is masked from 'package:MASS':
## 
##     select

## The following objects are masked from 'package:plyr':
## 
##     arrange, mutate, rename, summarise

## The following object is masked from 'package:stats':
## 
##     filter

## The following object is masked from 'package:graphics':
## 
##     layout

library(ltm)

## Loading required package: msm

## 
## Attaching package: 'msm'

## The following object is masked from 'package:boot':
## 
##     cav

## Loading required package: polycor

## 
## Attaching package: 'polycor'

## The following object is masked from 'package:psych':
## 
##     polyserial

## 
## Attaching package: 'ltm'

## The following object is masked from 'package:psych':
## 
##     factor.scores

library(polycor)
library(descr)
library(vcd)

## Loading required package: grid

library(leaps)
library(relaimpo)

## Loading required package: survey

## Loading required package: Matrix

## 
## Attaching package: 'survey'

## The following object is masked from 'package:graphics':
## 
##     dotchart

## Loading required package: mitools

## This is the global version of package relaimpo.

## If you are a non-US user, a version with the interesting additional metric pmvd is available

## from Ulrike Groempings web site at prof.beuth-hochschule.de/groemping.

library(FSelector)
library(caret)

## Warning: package 'caret' was built under R version 3.6.2

## 
## Attaching package: 'caret'

## The following object is masked from 'package:survival':
## 
##     cluster

library(factoextra)

## Warning: package 'factoextra' was built under R version 3.6.2

## Welcome! Want to learn more? See two factoextra-related books at https://goo.gl/ve3WBa

library(FactoMineR)

## Warning: package 'FactoMineR' was built under R version 3.6.2

library(nFactors)

## Warning: package 'nFactors' was built under R version 3.6.2

## 
## Attaching package: 'nFactors'

## The following object is masked from 'package:lattice':
## 
##     parallel

library(cluster)
library(stats)
library(forecast)

## Warning: package 'forecast' was built under R version 3.6.2

## Registered S3 method overwritten by 'xts':
##   method     from
##   as.zoo.xts zoo

## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo

## Registered S3 methods overwritten by 'forecast':
##   method             from    
##   fitted.fracdiff    fracdiff
##   residuals.fracdiff fracdiff

## 
## Attaching package: 'forecast'

## The following object is masked from 'package:ggpubr':
## 
##     gghistogram

## The following object is masked from 'package:rcompanion':
## 
##     accuracy

library(Metrics)

## Warning: package 'Metrics' was built under R version 3.6.2

## 
## Attaching package: 'Metrics'

## The following object is masked from 'package:forecast':
## 
##     accuracy

## The following objects are masked from 'package:caret':
## 
##     precision, recall

## The following object is masked from 'package:rcompanion':
## 
##     accuracy

## The following object is masked from 'package:FSA':
## 
##     se

library(tsoutliers)

## Warning: package 'tsoutliers' was built under R version 3.6.2

library(bfast)

## Warning: package 'bfast' was built under R version 3.6.2

Cargamos el fichero de datos y lo convertimos a dataframe.

dfMedidas <- read.csv("earthquake.csv", sep = ",", dec = ".", header = TRUE)

dfMedidas

Tenemos datos de las siguientes medidas de terremotos:

• md: Magnitud de duración.

• richter: Escala de richter.

• mw: Magnitud de energía.

• ms: Magnitud de ondas superficiales.

• mb: Magnitud de ondas de volumen.

• xm: max(md,richter,mw,ms,mb)

Vamos a usar siempre que sea posible la medida xm, ya que es el valor máximo de las demás medidas, y siempre hay datos para ella.

Calculamos la media, mediana, desviación estándar, valor máximo y mínimo, cuartiles y error estandar de la media para los valores de xm.

mean(dfMedidas$xm) #Cálculo de la Media para xm 

## [1] 4.056038

median(dfMedidas$xm) #Cálculo de la Mediana para xm

## [1] 3.9

sd(dfMedidas$xm) #Cálculo de la Desviación estándar para xm

## [1] 0.574085

range(dfMedidas$xm) #Valor máximo y mínimo de xm

## [1] 3.5 7.9

quantile(dfMedidas$xm) #Cuartiles para xm

##   0%  25%  50%  75% 100% 
##  3.5  3.6  3.9  4.4  7.9

sd(dfMedidas$xm) /sqrt(length(dfMedidas$xm)) #Error estándar de la media para mx

## [1] 0.003705162

Calculamos la Media de los valores de xm en cada país o región.

tapply(dfMedidas$xm, dfMedidas$country, mean)

##           #NAME?        aegeansea          albania       azerbaijan 
##         3.860000         3.913558         4.150000         4.613333 
##         blacksea         bulgaria     cyprus_greek   cyprus_turkish 
##         4.426667         4.450568         4.123308         3.942308 
##            egypt         georgia            greece             iran 
##         4.400000         4.472050         4.089719         4.578902 
##             iraq           israel        macedonia    mediterranean 
##         4.372131         3.800000         4.182143         4.127173 
##          romania           russia            syria           turkey 
##         4.400000         4.521122         3.956494         3.979249 
##  turkiye_armenia turkiye_georgia      turkiye_iran     turkiye_iraq 
##         4.620000         3.700000         4.112162         4.010000 
##    turkiye_syria         ukrainia 
##         3.933333         4.500000

Calculamos el número de terremotos que hay en cada país o región, junto a la media, desviación estándar, mediana, valor mínimo y máximo, y cuartiles del valor de la medida xm para cada país o región.

Summarize(xm ~ country, data = dfMedidas)

Calculamos esas medidas y alguna más pero con describeBy.

describeBy(dfMedidas$xm, group = dfMedidas$country, digits= 4)

## 
##  Descriptive statistics by group 
## group: #NAME?
##    vars  n mean   sd median trimmed  mad min max range skew kurtosis   se
## X1    1 10 3.86 0.28    3.9    3.85 0.37 3.5 4.3   0.8 0.07    -1.52 0.09
## -------------------------------------------------------- 
## group: aegeansea
##    vars    n mean  sd median trimmed mad min max range skew kurtosis   se
## X1    1 1748 3.91 0.5    3.7    3.82 0.3 3.5 7.2   3.7 2.11     6.23 0.01
## -------------------------------------------------------- 
## group: albania
##    vars n mean   sd median trimmed  mad min max range skew kurtosis   se
## X1    1 2 4.15 0.64   4.15    4.15 0.67 3.7 4.6   0.9    0    -2.75 0.45
## -------------------------------------------------------- 
## group: azerbaijan
##    vars   n mean  sd median trimmed  mad min max range skew kurtosis   se
## X1    1 150 4.61 0.6   4.65    4.59 0.52 3.5 6.7   3.2 0.42     0.26 0.05
## -------------------------------------------------------- 
## group: blacksea
##    vars  n mean   sd median trimmed  mad min max range skew kurtosis   se
## X1    1 90 4.43 0.74    4.3    4.35 0.74 3.5 6.8   3.3 0.79     0.02 0.08
## -------------------------------------------------------- 
## group: bulgaria
##    vars   n mean   sd median trimmed  mad min max range skew kurtosis   se
## X1    1 176 4.45 0.72    4.5     4.4 0.89 3.5 7.1   3.6 0.81     0.91 0.05
## -------------------------------------------------------- 
## group: cyprus_greek
##    vars   n mean   sd median trimmed  mad min max range skew kurtosis   se
## X1    1 133 4.12 0.57      4    4.04 0.44 3.5 6.3   2.8 1.35     1.52 0.05
## -------------------------------------------------------- 
## group: cyprus_turkish
##    vars  n mean   sd median trimmed  mad min max range skew kurtosis   se
## X1    1 26 3.94 0.61    3.7    3.83 0.22 3.5 6.5     3 2.85     8.89 0.12
## -------------------------------------------------------- 
## group: egypt
##    vars n mean   sd median trimmed  mad min max range skew kurtosis  se
## X1    1 2  4.4 0.42    4.4     4.4 0.44 4.1 4.7   0.6    0    -2.75 0.3
## -------------------------------------------------------- 
## group: georgia 
##    vars   n mean   sd median trimmed  mad min max range skew kurtosis   se
## X1    1 322 4.47 0.63    4.5    4.44 0.59 3.5 6.9   3.4 0.62     0.56 0.03
## -------------------------------------------------------- 
## group: greece
##    vars    n mean   sd median trimmed  mad min max range skew kurtosis
## X1    1 3560 4.09 0.57    3.9    4.01 0.44 3.5 7.7   4.2 1.36      2.2
##      se
## X1 0.01
## -------------------------------------------------------- 
## group: iran
##    vars   n mean   sd median trimmed  mad min max range skew kurtosis   se
## X1    1 346 4.58 0.54    4.6    4.57 0.44 3.5 7.1   3.6 0.53      1.5 0.03
## -------------------------------------------------------- 
## group: iraq
##    vars   n mean   sd median trimmed  mad min max range skew kurtosis   se
## X1    1 122 4.37 0.59   4.35    4.34 0.67 3.5 5.9   2.4 0.34    -0.56 0.05
## -------------------------------------------------------- 
## group: israel
##    vars n mean sd median trimmed mad min max range skew kurtosis se
## X1    1 1  3.8 NA    3.8     3.8   0 3.8 3.8     0   NA       NA NA
## -------------------------------------------------------- 
## group: macedonia
##    vars  n mean   sd median trimmed  mad min max range skew kurtosis  se
## X1    1 28 4.18 0.51      4    4.16 0.52 3.5 5.3   1.8 0.51    -1.04 0.1
## -------------------------------------------------------- 
## group: mediterranean
##    vars    n mean   sd median trimmed  mad min max range skew kurtosis
## X1    1 4843 4.13 0.57      4    4.05 0.59 3.5 7.4   3.9 1.27     1.85
##      se
## X1 0.01
## -------------------------------------------------------- 
## group: romania
##    vars  n mean   sd median trimmed  mad min max range skew kurtosis  se
## X1    1 44  4.4 0.69    4.4    4.34 0.74 3.5   6   2.5 0.66     -0.5 0.1
## -------------------------------------------------------- 
## group: russia
##    vars   n mean  sd median trimmed  mad min max range skew kurtosis   se
## X1    1 303 4.52 0.5    4.5    4.49 0.44 3.6 6.5   2.9 0.69     0.84 0.03
## -------------------------------------------------------- 
## group: syria
##    vars   n mean   sd median trimmed mad min max range skew kurtosis   se
## X1    1 154 3.96 0.45    3.8    3.89 0.3 3.5 5.8   2.3 1.39      1.9 0.04
## -------------------------------------------------------- 
## group: turkey
##    vars     n mean   sd median trimmed  mad min max range skew kurtosis
## X1    1 11850 3.98 0.55    3.8    3.88 0.44 3.5 7.9   4.4 1.69     3.47
##      se
## X1 0.01
## -------------------------------------------------------- 
## group: turkiye_armenia
##    vars n mean   sd median trimmed mad min max range  skew kurtosis   se
## X1    1 5 4.62 0.28    4.7    4.62 0.3 4.2 4.9   0.7 -0.44     -1.7 0.12
## -------------------------------------------------------- 
## group: turkiye_georgia 
##    vars n mean sd median trimmed mad min max range skew kurtosis se
## X1    1 1  3.7 NA    3.7     3.7   0 3.7 3.7     0   NA       NA NA
## -------------------------------------------------------- 
## group: turkiye_iran
##    vars  n mean   sd median trimmed  mad min max range skew kurtosis   se
## X1    1 74 4.11 0.65    4.1    4.02 0.59 3.5 7.6   4.1 2.34     9.32 0.08
## -------------------------------------------------------- 
## group: turkiye_iraq
##    vars  n mean   sd median trimmed mad min max range skew kurtosis   se
## X1    1 10 4.01 0.61   3.75     3.9 0.3 3.5 5.4   1.9 1.16    -0.02 0.19
## -------------------------------------------------------- 
## group: turkiye_syria
##    vars n mean   sd median trimmed  mad min max range skew kurtosis   se
## X1    1 6 3.93 0.41   3.85    3.93 0.22 3.5 4.7   1.2 0.84     -0.8 0.17
## -------------------------------------------------------- 
## group: ukrainia
##    vars n mean sd median trimmed mad min max range skew kurtosis se
## X1    1 1  4.5 NA    4.5     4.5   0 4.5 4.5     0   NA       NA NA

Calculamos la media de los valores de xm con un intervalo de confianza del 95%

CI(dfMedidas$xm, ci = 0.95)

##    upper     mean    lower 
## 4.063300 4.056038 4.048775

calculamos la desviación absoluta media, la desviación estándar y el rango intercuartil para los valores de xm, respecto a albania, grecia, mediterraneo y turquía.

mad(dfMedidas[which(dfMedidas$country == "albania"), "xm"]) #desviación absoluta media respecto a albania

## [1] 0.66717

mad(dfMedidas[which(dfMedidas$country == "greece"), "xm"]) #desviación absoluta media respecto a grecia

## [1] 0.44478

mad(dfMedidas[which(dfMedidas$country == "mediterranean"), "xm"]) #desviación absoluta media respecto al mediterraneo

## [1] 0.59304

mad(dfMedidas[which(dfMedidas$country == "turkey"), "xm"]) #desviación absoluta media respecto a turquia

## [1] 0.44478

sd(dfMedidas[which(dfMedidas$country == "albania"), "xm"]) #desviación estándar respecto a albania

## [1] 0.6363961

sd(dfMedidas[which(dfMedidas$country == "greece"), "xm"]) #desviación estándar respecto a grecia

## [1] 0.5689156

sd(dfMedidas[which(dfMedidas$country == "mediterranean"), "xm"]) #desviación estándar respecto al mediterraneo

## [1] 0.5738599

sd(dfMedidas[which(dfMedidas$country == "turkey"), "xm"]) #desviación estándar respecto a turquia

## [1] 0.5507318

IQR(dfMedidas[which(dfMedidas$country == "albania"), "xm"]) #rango intercuartil respecto a albania

## [1] 0.45

IQR(dfMedidas[which(dfMedidas$country == "greece"), "xm"]) #rango intercuartil respecto a grecia

## [1] 0.8

IQR(dfMedidas[which(dfMedidas$country == "mediterranean"), "xm"]) #rango intercuartil respecto al mediterraneo

## [1] 0.7

IQR(dfMedidas[which(dfMedidas$country == "turkey"), "xm"]) #rango intercuartil respecto a turquia

## [1] 0.6

Comprobamos que tal se ajustan los valores de xm a una función normal para albania, grecia, mediterraneo y turquía.

plotNormalHistogram(dfMedidas[which(dfMedidas$country == "albania"), "xm"])

plotNormalHistogram(dfMedidas[which(dfMedidas$country == "greece"), "xm"])

plotNormalHistogram(dfMedidas[which(dfMedidas$country == "mediterranean"), "xm"])

plotNormalHistogram(dfMedidas[which(dfMedidas$country == "turkey"), "xm"])

Como se puede ver no se pueden ajustar a una normal.

Intentamos conseguir una mayor aproximación a los datos de xm para el mediterráneo, por medio de asimetría estadística (skewness).

op <- par(mar = c(3, 3, 4, 2) + 0.1)

hist(dfMedidas[which(dfMedidas$country == "mediterranean"), "xm"], col = "light blue", probability=TRUE, main = paste("skewness =", round(skewness(dfMedidas[which(dfMedidas$country == "mediterranean"), "xm"]), digits=2)),xlab="", ylab="")
lines(density(dfMedidas[which(dfMedidas$country == "mediterranean"), "xm"]), col ="red", lwd = 5)
rug(dfMedidas[which(dfMedidas$country == "mediterranean"), "xm"], col ="green", lwd = 3)

par(op)

skewness(dfMedidas[which(dfMedidas$country == "mediterranean"), "xm"])

## [1] 1.27069

Se aprecia una mejor aproximación.

Por medio de la medida de curtosis vemos si los datos se extienden más o menos que la normal en las colas y si hay mayor densidad alrededor del valor medio. En este caso seguimos viéndolo para los terremotos del mediterraneo.

xm_mediterranean <- dfMedidas[which(dfMedidas$country == "mediterranean"), "xm"]

qqnorm(xm_mediterranean, main = paste("curtosis =", round(kurtosis(xm_mediterranean), digits = 2), "(gaussian)"))
qqline(xm_mediterranean, col = "red")
op <- par(fig=c(0.02, 0.5, 0.5, 0.98), new = TRUE)
hist(xm_mediterranean, probability=T, col="light blue", xlab = "", ylab = "", main = "", axes = FALSE)
lines(density(xm_mediterranean), col = "red", lwd=2)
box()

par(op)

Comprobamos a que función de distribución se asemejan más nuestros datos de los terremotos con el valor de xm ocurridos en el mediterraneo.

descdist(xm_mediterranean, discrete = FALSE)

## summary statistics
## ------
## min:  3.5   max:  7.4 
## median:  4 
## mean:  4.127173 
## estimated sd:  0.5738599 
## estimated skewness:  1.271477 
## estimated kurtosis:  4.850741

Como se puede apreciar se acerca bastante a la función gamma, por lo que se podría aproximar a dicha función, o bien a la función beta.

Comprobamos la bondad del ajuste a una función normal de nuestros datos en el mediterráneo, mediante el test de Shapiro-Wilk.

shapiro.test(xm_mediterranean)

## 
##  Shapiro-Wilk normality test
## 
## data:  xm_mediterranean
## W = 0.88511, p-value < 2.2e-16

Como se puede apreciar el p-value está muy por debajo de 0,05, por lo que no se puede considerar que se asemeje en absoluto a una normal.

Hacemos una gráfica Q-Q de los datos del mediterráneo para ver visualmente si se puede asemejar a una normal, y vemos definitivamente que no, ya que los valores casi nunca entran en el intervalo de confianza.

ggqqplot(xm_mediterranean) #gráfica Q-Q

Hacemos otra gráfica Q-Q, pero en este caso para todas las regiones o países.

qqPlot(dfMedidas$xm) #La misma gráfica Q-Q pero con otra librería y en general 

## [1]  6718 10319

Se observa que los valores tampoco se pueden asemejar a una normal.

Para más seguridad, realizamos las pruebas de Anderson-Darling, Cramér-von Mises y Kolmogorov- Smirnov para los terremotos acontecidos en el mediterráneo por un lado, y en todos los países o regiones por otro.

ad.test(xm_mediterranean)

## 
##  Anderson-Darling normality test
## 
## data:  xm_mediterranean
## A = 147.06, p-value < 2.2e-16

cvm.test(xm_mediterranean)

## Warning in cvm.test(xm_mediterranean): p-value is smaller than 7.37e-10,
## cannot be computed more accurately

## 
##  Cramer-von Mises normality test
## 
## data:  xm_mediterranean
## W = 23.801, p-value = 7.37e-10

lillie.test(xm_mediterranean)

## 
##  Lilliefors (Kolmogorov-Smirnov) normality test
## 
## data:  xm_mediterranean
## D = 0.13818, p-value < 2.2e-16

ad.test(dfMedidas$xm)

## 
##  Anderson-Darling normality test
## 
## data:  dfMedidas$xm
## A = 1010, p-value < 2.2e-16

cvm.test(dfMedidas$xm)

## Warning in cvm.test(dfMedidas$xm): p-value is smaller than 7.37e-10, cannot
## be computed more accurately

## 
##  Cramer-von Mises normality test
## 
## data:  dfMedidas$xm
## W = 171.41, p-value = 7.37e-10

lillie.test(dfMedidas$xm)

## 
##  Lilliefors (Kolmogorov-Smirnov) normality test
## 
## data:  dfMedidas$xm
## D = 0.16638, p-value < 2.2e-16

Todas las pruebas indican que los valores no se pueden aproximar a una normal, al tener un p-value mucho menor de 0,05.

Determinamos tres tipos de alertas (amarilla, naranja y roja) y discretizamos los valores de xm en el mediterráneo.

xm_mediterranean_alert <- classIntervals(xm_mediterranean, n = 3, style = "fixed",fixedBreaks = c(2.9, 4.5, 6.0, 7.9), intervalClosure = "left") # Frecuencias de las alertas
plot(xm_mediterranean_alert, pal = c("yellow", "orange", "red"), main = "xm_mediterranean")

Vemos una posibilidad de 151,532,656,696 posibles particiónes de los valores de xm en intervalos.

(classIntervals(dfMedidas$xm, style = "fisher", intervalClosure = "left"))

## style: fisher
##   one of 151,532,656,696 possible partitions of this variable into 16 classes
##  [3.5,3.55) [3.55,3.65) [3.65,3.75) [3.75,3.85) [3.85,3.95) [3.95,4.05) 
##        4071        3196        2506        1992        1559        1403 
## [4.05,4.25) [4.25,4.45) [4.45,4.65) [4.65,4.85) [4.85,5.05) [5.05,5.35) 
##        2142        1907        1414        1295         934         765 
## [5.35,5.65) [5.65,6.05) [6.05,6.65)  [6.65,7.9] 
##         454         222          90          57

Vemos Una de las 903 posibles particiones de los datos de mx en 3 intervalos.

(classIntervals(dfMedidas$xm, n = 3, style = "quantile", intervalClosure = "left"))

## style: quantile
##   one of 903 possible partitions of this variable into 3 classes
## [3.5,3.7) [3.7,4.2) [4.2,7.9] 
##      7267      8644      8096

Dividimos la profundidad en 5 intervalos después de hacer varias pruebas y determinar cómo óptimo dicho número de particiones, y añadimos una nueva columna con esos intervalos.

depthInterval <- dfMedidas$depth

depthInterval[dfMedidas$depth < 7] <- "[0,7)"
depthInterval[dfMedidas$depth >= 7 & dfMedidas$depth < 15] <- "[7,15)"
depthInterval[dfMedidas$depth >= 15 & dfMedidas$depth < 36] <- "[15,36)"
depthInterval[dfMedidas$depth >= 36 & dfMedidas$depth < 78] <- "[36,78)"
depthInterval[dfMedidas$depth >= 78 & dfMedidas$depth <= 225] <- "[78,225]"

dfMedidas <- cbind(dfMedidas[1:10],depthInterval,dfMedidas[11:17])

Hacemos un histograma para cada intervalo de profundidad, para ver que tal se aproximan los valores de la medida xm a una normal teniendo en cuenta su profundidad.

plotNormalHistogram(dfMedidas[which(dfMedidas$depthInterval == "[0,7)"), "xm"])

plotNormalHistogram(dfMedidas[which(dfMedidas$depthInterval == "[7,15)"), "xm"])

plotNormalHistogram(dfMedidas[which(dfMedidas$depthInterval == "[15,36)"), "xm"])

plotNormalHistogram(dfMedidas[which(dfMedidas$depthInterval == "[36,78)"), "xm"])

plotNormalHistogram(dfMedidas[which(dfMedidas$depthInterval == "[78,225]"), "xm"])

Se puede apreciar que los datos no se asemejan nada o casi nada a una normal en ningún intervalo.

Vemos a que función o funciones se pueden aproximar los datos de cada intervalo:

Se aprecia que en el intervalo de profundidad entre 0 y 7, los valores se pueden aproximar muy bien a una función gamma.

descdist(dfMedidas[which(dfMedidas$depthInterval == "[0,7)"), "xm"])

## summary statistics
## ------
## min:  3.5   max:  7.2 
## median:  3.7 
## mean:  3.845054 
## estimated sd:  0.4064809 
## estimated skewness:  2.130021 
## estimated kurtosis:  9.871725

Se aprecia que en el intervalo de profundidad entre 7 y 15, los valores se pueden aproximar a una función beta y esta algo cerca de una gamma.

descdist(dfMedidas[which(dfMedidas$depthInterval == "[7,15)"), "xm"])

## summary statistics
## ------
## min:  3.5   max:  7.5 
## median:  3.8 
## mean:  3.959905 
## estimated sd:  0.5163223 
## estimated skewness:  1.917332 
## estimated kurtosis:  8.17031

Se aprecia que en el intervalo de profundidad entre 15 y 36, los valores también se pueden aproximar una función beta y algo a una gamma y está más cerca que el anterior de una lognormal.

descdist(dfMedidas[which(dfMedidas$depthInterval == "[15,36)"), "xm"])

## summary statistics
## ------
## min:  3.5   max:  7.9 
## median:  4.2 
## mean:  4.270902 
## estimated sd:  0.6432241 
## estimated skewness:  0.881727 
## estimated kurtosis:  3.756313

Se aprecia que en el intervalo de profundidad entre 36 y 78, los valores se pueden aproximar una función gamma y a una lognormal, y es el que está más cerca de una normal.

descdist(dfMedidas[which(dfMedidas$depthInterval == "[36,78)"), "xm"])

## summary statistics
## ------
## min:  3.5   max:  7.6 
## median:  4.5 
## mean:  4.496967 
## estimated sd:  0.5934879 
## estimated skewness:  0.5523241 
## estimated kurtosis:  3.553707

Se aprecia que en el intervalo de profundidad entre 78 y 225, los valores se pueden aproximar una función beta y se acercan a la gamma.

descdist(dfMedidas[which(dfMedidas$depthInterval == "[78,225]"), "xm"])

## summary statistics
## ------
## min:  3.5   max:  7.7 
## median:  4.3 
## mean:  4.371499 
## estimated sd:  0.6736299 
## estimated skewness:  0.9394736 
## estimated kurtosis:  4.051817

Hacemos un diagrama de dispersión con diferentes técnicas, de xm en función de la profundidad.

par(mfrow=c(1,1))
plot(dfMedidas$depth, dfMedidas$xm, main="Diagrama de Dispersión", xlab="Profundidad en metros", ylab="xm", pch=19)
abline(lm(dfMedidas$xm ~ dfMedidas$depth), col = "red") # regresión (y~x)

smoothScatter(dfMedidas$depth, dfMedidas$xm,xlab="Profundidad en metros", ylab="xm", cex= 2,colramp = colorRampPalette(c("white", "blue"), space = "Lab"))

scatterplot(xm ~ depth, data = dfMedidas, main ="Diagrama de Dispersión", xlab ="Profundidad en metros", ylab ="xm", ellipse = TRUE, lwd = 2, lty= 2)

scatter.hist(dfMedidas$depth, dfMedidas$xm, main ="Diagrama de Dispersión", xlab ="Profundidad en metros", ylab ="xm")

#Igual que la anterior, pero con histogramas y coeficiente de regresión

Se puede apreciar que por lo general, la mayoría de terremotos se producen a poca o media profundidad, y en la mayoría de los casos tienen una medida xm menor de ~ 5.5.

Correlación entre todas las medidas de terremotos (Sin tener en cuenta los ceros).

scatterplotMatrix(~ md[which((dfMedidas$md>0)&(dfMedidas$richter>0)&(dfMedidas$mw>0)&(dfMedidas$ms>0)&(dfMedidas$mb>0))] + richter[which((dfMedidas$md>0)&(dfMedidas$richter>0)&(dfMedidas$mw>0)&(dfMedidas$ms>0)&(dfMedidas$mb>0))] + mw[which((dfMedidas$md>0)&(dfMedidas$richter>0)&(dfMedidas$mw>0)&(dfMedidas$ms>0)&(dfMedidas$mb>0))] + ms[which((dfMedidas$md>0)&(dfMedidas$richter>0)&(dfMedidas$mw>0)&(dfMedidas$ms>0)&(dfMedidas$mb>0))] + mb[which((dfMedidas$md>0)&(dfMedidas$richter>0)&(dfMedidas$mw>0)&(dfMedidas$ms>0)&(dfMedidas$mb>0))], data = dfMedidas)

Como se puede apreciar hay bastante correlación entre todas ellas, lo que nos reafirma en nuestra decisión de usar siempre la medida xm.

Calculamos la covarianza y correlación de la profundidad con la medida xm.

cov(dfMedidas$depth, dfMedidas$xm)

## [1] 4.037828

cor(dfMedidas$depth, dfMedidas$xm)

## [1] 0.3029259

Se aprecia poca correlación, ya que esta es menor de 0.4, que es cuando se empezaría a considerar moderada, aunque se puede apreciar que cuando la profundidad aumenta, la escala xm de los terremotos que se producen a esa profundidad, también aumenta al ser positiva la covarianza.

Se usa la función hetcor para calcular las correlaciones de xm con la latitud, longitud y profundidad, usando para cada una el tipo de correlación que le corresponde, siendo usada en este caso siempre la de Pearson.

hetcor(dfMedidas$xm, dfMedidas$lat, dfMedidas$long, dfMedidas$depth, dfMedidas$dist)

## 
## Two-Step Estimates
## 
## Correlations/Type of Correlation:
##                 dfMedidas$xm dfMedidas.lat dfMedidas.long dfMedidas.depth
## dfMedidas$xm               1       Pearson        Pearson         Pearson
## dfMedidas.lat         0.0747             1        Pearson         Pearson
## dfMedidas.long       0.06665       0.07588              1         Pearson
## dfMedidas.depth       0.3725      -0.09832     -0.0003035               1
## dfMedidas.dist      0.004511       0.07934        0.02833         0.02869
##                 dfMedidas.dist
## dfMedidas$xm           Pearson
## dfMedidas.lat          Pearson
## dfMedidas.long         Pearson
## dfMedidas.depth        Pearson
## dfMedidas.dist               1
## 
## Standard Errors:
##                 dfMedidas$xm dfMedidas.lat dfMedidas.long dfMedidas.depth
## dfMedidas$xm                                                             
## dfMedidas.lat       0.009914                                             
## dfMedidas.long      0.009925      0.009912                               
## dfMedidas.depth     0.008586      0.009873        0.00997                
## dfMedidas.dist      0.009969      0.009907       0.009962        0.009961
## 
## n = 10062 
## 
## P-values for Tests of Bivariate Normality:
##                 dfMedidas$xm dfMedidas.lat dfMedidas.long dfMedidas.depth
## dfMedidas$xm                                                             
## dfMedidas.lat     3.953e-210                                             
## dfMedidas.long             0             0                               
## dfMedidas.depth            0             0              0                
## dfMedidas.dist             0             0              0               0

Se puede apreciar que la correlación con las demás variables es incluso peor que con la profundidad.

Creamos diferentes intervalos de xm en función de lo que pensamos que consideramos como un terremoto débil, uno medio o uno fuerte, y añadimos una columna en el dataframe con dichos intervalos.

xm_Interval <- dfMedidas$xm

xm_Interval[dfMedidas$xm < 4.5] <- "[3.5,4.5)"
xm_Interval[dfMedidas$xm >= 4.5 & dfMedidas$xm < 6] <- "[4.5,6)"
xm_Interval[dfMedidas$xm >= 6 & dfMedidas$xm <= 7.9] <- "[6,7.9]"

dfMedidas <- cbind(dfMedidas[1:12],xm_Interval,dfMedidas[13:18])

Vemos la proporción de terremotos que hay en cada intervalo de xm, y la proporción de cada intervalo de xm en función de cada intervalo de profundidad.

freq(xm_Interval, plot=FALSE) #frecuencias marginales

## xm_Interval 
##           Frequency  Percent
## [3.5,4.5)     18776  78.2105
## [4.5,6)        5044  21.0105
## [6,7.9]         187   0.7789
## Total         24007 100.0000

round(prop.table(table(depthInterval, xm_Interval), 2), 2) #tabla de contingencia

##              xm_Interval
## depthInterval [3.5,4.5) [4.5,6) [6,7.9]
##      [0,7)         0.38    0.13    0.11
##      [15,36)       0.17    0.34    0.33
##      [36,78)       0.06    0.24    0.20
##      [7,15)        0.36    0.23    0.28
##      [78,225]      0.03    0.06    0.09

Obtenemos dos tablas, una con el número de terremotos esperados por chi cuadrado en cada intervalo y otra con el numero real de terremotos

round(chisq.test(table(depthInterval, xm_Interval))$expected, 0) #esperada

##              xm_Interval
## depthInterval [3.5,4.5) [4.5,6) [6,7.9]
##      [0,7)         6057    1627      60
##      [15,36)       3865    1038      38
##      [36,78)       1934     520      19
##      [7,15)        6289    1689      63
##      [78,225]       631     170       6

table(depthInterval, xm_Interval) #real

##              xm_Interval
## depthInterval [3.5,4.5) [4.5,6) [6,7.9]
##      [0,7)         7081     643      20
##      [15,36)       3146    1734      62
##      [36,78)       1216    1220      37
##      [7,15)        6846    1143      52
##      [78,225]       487     304      16

Hacemos el test de chi cuadrado y luego representamos los valores que nos da, viendo si los valores aciertan o fallan respecto a lo esperado.

chisq.test(table(depthInterval, xm_Interval))

## 
##  Pearson's Chi-squared test
## 
## data:  table(depthInterval, xm_Interval)
## X-squared = 3019.2, df = 8, p-value < 2.2e-16

El valor de p-value es < 2.2e-16, por lo que es mucho menor de 0.05 y por tanto se puede rechazar la hipótesis nula, y afirmar que existe una relación significativa entre xm y la profundidad.

Obtenemos una tabla en la que se puede obsevar la probabilidad de que un terremoto tenga una escala perteneciente a un intervalo determinado, sabiendo que se produce en un determinado intervalo de profundidad.

assoc(table(depthInterval, xm_Interval), shade=TRUE)

Se observa que:

• Para profundidad entre 0 y 7, es bastante probable que se produzcan terremotos de entre 3.5 y 4.5, muy probale que no se produzcan de medidas entre 4.5 y 6, y probablemente tampoco se producirán de entre 6 y 7.9.

• Para profundidad entre 7 y 15, tenemos un escenario parecido al anterior.

• Para profundidad entre 15 y 36, es bastante probable que no se produzcan terremotos de entre 3.5 y 4.5, muy probale que se produzcan de medidas entre 4.5 y 6, y algo probable que se produzcan algunos de entre 6 y 7.9.

• Para profundidad entre 36 y 78, es bastante probable que no se produzcan terremotos de entre 3.5 y 4.5, muy probale que se produzcan de medidas entre 4.5 y 6, y es en el intervalo de produndidad donde es más probable que se produzcan terremotos de entre 6 y 7.9, aunque sigue siendo poca, al haber tan poca cantidad de de terremotos con esas medidas.

• Para profundidad entre 78 y 225, es probable que no se produzcan terremotos de entre 3.5 y 4.5, bastante probable que se produzcan de medidas entre 4.5 y 6, y algo probable que se produzcan algunos de entre 6 y 7.9.

Creamos 2 nuevas columnas: una para guardar los años, y otra para guardar los meses.

year <- substring(dfMedidas$date, 1, 4)

dfMedidas <- cbind(dfMedidas[1],year,dfMedidas[2:19])

month <- substring(dfMedidas$date, 6, 7)

dfMedidas <- cbind(dfMedidas[1:2],month,dfMedidas[3:20])

Dibujamos un histograma y una distribución de densidad con los diferentes valores de xm, poniendo la línea roja en 6, que es lo que consideramos un terremoto grave.

qplot(dfMedidas$xm, data=dfMedidas, geom="density", alpha=I(.5)) + xlab(bquote(xm)) + geom_histogram(fill = "#4271AE", colour = "#1F3552", alpha = 0.6) + scale_x_continuous(breaks = seq(0, 10, 0.5), limits=c(2, 8)) + geom_vline(xintercept=6, linetype="dashed", color = "red") + ggtitle("Histograma Escala de xm")

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## Warning: Removed 2 rows containing missing values (geom_bar).

qplot(xm, data = dfMedidas, geom = "density", alpha = I(.5)) + labs(x = bquote(xm)) + geom_density(fill = "#4271AE", colour = "#1F3552", alpha = 0.6) + scale_x_continuous(breaks = seq(0, 10, 0.5), limits = c(2, 8)) + geom_vline(xintercept = 6, linetype="dashed", color = "red") + ggtitle("Distribución de Densidad Escala de xm")

Se aprecia que a medida que la escala del terremoto aumenta el número de terremotos de dicha escala disminuye, y viceversa.

Se estandarizan las diferentes medidas de teremotos (excepto mw, ya que tiene muchos Na) y se comparan sus histogramas.

h1 <- ggplot(data=dfMedidas, aes((xm - mean(xm)) / sd(xm))) + geom_histogram(breaks=seq(-3, 3, by = 0.5), col = "gray", fill = "green", alpha = 0.2) + labs(title = "Histograma") + labs(x = bquote(xm), y = "Cantidad")

h2 <- ggplot(data=dfMedidas, aes((md - mean(md)) / sd(md))) + geom_histogram(breaks=seq(-3, 3, by = 0.5), col = "gray", fill = "green", alpha = 0.2) + labs(title = "Histograma") + labs(x = bquote("md (Magnitud de duración)"), y = "Cantidad")

h3 <- ggplot(data=dfMedidas, aes((richter - mean(richter)) / sd(richter))) + geom_histogram(breaks=seq(-3, 3, by = 0.5), col = "gray", fill = "green", alpha = 0.2) + labs(title = "Histograma") + labs(x = bquote(richter), y = "Cantidad")

h4 <- ggplot(data=dfMedidas, aes((ms - mean(ms)) / sd(ms))) + geom_histogram(breaks=seq(-3, 3, by = 0.5), col = "gray", fill = "green", alpha = 0.2) + labs(title = "Histograma") + labs(x = bquote("ms (Magnitud de ondas superficiales)"), y = "Cantidad")

h5 <- ggplot(data=dfMedidas, aes((mb - mean(mb)) / sd(mb))) + geom_histogram(breaks=seq(-3, 3, by = 0.5), col = "gray", fill = "green", alpha = 0.2) + labs(title = "Histograma") + labs(x = bquote("mb (Magnitud de ondas de volumen)"), y = "Cantidad")

grid.arrange(h1, h2, h3, h4, h5, ncol = 2, nrow = 3)

Ponemos las diferentes medidas de xm en latitud y longitud para hacer un mapa de puntos con terremotos, usando diferentes colores para valores de xm mayores y menores que 5.5, además de aumentar el tamaño de los puntos a medida que aumenta la profundidad.

ggplot(dfMedidas, aes(x = long, y = lat, size=depth, fill = xm < 5.5)) + ylab("Latitud") + xlab("Longitud") + geom_point(shape = 21) + ggtitle("Medida xm") + scale_size(range = c(1, 5))

Comparamos las diferentes medidas de fuerza de un terremoto con la profundidad.

ggplot(dfMedidas, aes(x = depth, y = xm, fill = xm < 5.5)) + ylab("xm") + xlab("Profundidad") + geom_point(shape = 21) + ggtitle("Medida xm v Profundidad") + geom_hline(yintercept=6, linetype="dashed", color = "red")

ggplot(dfMedidas, aes(x = depth, y = md, fill = md < 5.5)) + ylab("md") + xlab("Profundidad") + geom_point(shape = 21) + ggtitle("Magnitud de duración (md) v Profundidad") + geom_hline(yintercept=6, linetype="dashed", color = "red")

ggplot(dfMedidas, aes(x = depth, y = richter, fill = richter < 5.5)) + ylab("Escala de richter") + xlab("Profundidad") + geom_point(shape = 21) + ggtitle("Escala de richter v Profundidad") + geom_hline(yintercept=6, linetype="dashed", color = "red")

ggplot(dfMedidas, aes(x = depth, y = ms, fill = ms < 5.5)) + ylab("ms") + xlab("Profundidad") + geom_point(shape = 21) + ggtitle("Magnitud de ondas superficiales (ms) v Profundidad") + geom_hline(yintercept=6, linetype="dashed", color = "red")

ggplot(dfMedidas, aes(x = depth, y = mb, fill = mb < 5.5)) + ylab("mb") + xlab("Profundidad") + geom_point(shape = 21) + ggtitle("Magnitud de ondas de volumen (mb) v Profundidad") + geom_hline(yintercept=6, linetype="dashed", color = "red")

Como se ha visto anteriormente, se puede apreciar que por lo general, la mayoría de terremotos se producen a poca o media profundidad y en la mayoría de los casos tienen una medida xm menor de ~ 5.5

Pasamos “year” y “month” a numérico para poder compararlos.

dfMedidas$year <- as.numeric(year)
dfMedidas$month <- as.numeric(month)

Sacamos la medida de la importancia que tienen diferentes variables con respecto a xm, donde vemos que year es la más importante, pero cómo no debería tener mucho sentido saber el año en lo referente a predecir cuando se producirá un terremoto o de que magnitud será, consideraremos que la profundidad es la más importante.

importancia <- calc.relimp(lm(xm ~ year + month + depth + lat + long + dist, data =
dfMedidas), type = c("lmg", "last", "first", "betasq", "pratt", "genizi", "car"), rela = TRUE)

## Warning in rev(variances[[p]]) - variances[[p + 1]]: Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.

importancia

## Response variable: xm 
## Total response variance: 0.3145435 
## Analysis based on 10062 observations 
## 
## 6 Regressors: 
## year month depth lat long dist 
## Proportion of variance explained by model: 36.42%
## Metrics are normalized to sum to 100% (rela=TRUE). 
## 
## Relative importance metrics: 
## 
##                lmg         last        first       betasq         pratt
## year  0.6975144880 0.7688173336 6.645368e-01 0.7858954858  7.614332e-01
## month 0.0004082634 0.0005377904 9.093363e-05 0.0004427669  2.114169e-04
## depth 0.2373404122 0.1113401243 3.127421e-01 0.1087733714  1.943320e-01
## lat   0.0146015671 0.0113897045 1.257473e-02 0.0096766267  1.162255e-02
## long  0.0500043939 0.1076979223 1.000965e-02 0.0950322381  3.249638e-02
## dist  0.0001308755 0.0002171249 4.586455e-05 0.0001795111 -9.560359e-05
##             genizi          car
## year  0.6985512114 7.271511e-01
## month 0.0004234019 2.984050e-04
## depth 0.2368042635 2.190259e-01
## lat   0.0144045801 1.312157e-02
## long  0.0497078962 4.038965e-02
## dist  0.0001086470 1.338489e-05
## 
## Average coefficients for different model sizes: 
## 
##                  1X           2Xs           3Xs           4Xs
## year  -0.0159637418 -0.0157863597 -0.0156025431 -0.0154121667
## month  0.0010554052  0.0016052890  0.0019563902  0.0021279379
## depth  0.0136408031  0.0124558089  0.0111999339  0.0098680841
## lat    0.0338613864  0.0346111189  0.0341995879  0.0325998745
## long   0.0064482283  0.0089190976  0.0112270669  0.0133732894
## dist   0.0005365677  0.0001651884 -0.0002040547 -0.0005338381
##                 5Xs           6Xs
## year  -0.0152151498 -0.0150114553
## month  0.0021399327  0.0020137695
## depth  0.0084551858  0.0069562200
## lat    0.0297764687  0.0256851785
## long   0.0153579731  0.0171803321
## dist  -0.0007853369 -0.0009179032

plot(importancia)

Hacemos lo mismo que la importancia anterior pero basados en la profundidad, donde sale que mb es la mas importante para determinar la profundidad.

importancia <- calc.relimp(lm(depth ~ md + richter + mw + ms + mb, data =
dfMedidas), type = c("lmg", "last", "first", "betasq", "pratt", "genizi", "car"), rela = TRUE)

## Warning in rev(variances[[p]]) - variances[[p + 1]]: Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.

importancia

## Response variable: depth 
## Total response variance: 734.1677 
## Analysis based on 5003 observations 
## 
## 5 Regressors: 
## md richter mw ms mb 
## Proportion of variance explained by model: 10.43%
## Metrics are normalized to sum to 100% (rela=TRUE). 
## 
## Relative importance metrics: 
## 
##                lmg        last     first      betasq       pratt
## md      0.28179380 0.147450456 0.2758824 0.139358053  0.28957831
## richter 0.02696925 0.076254382 0.0360077 0.004799141 -0.01941411
## mw      0.10858376 0.133478189 0.1282842 0.011386866  0.05644515
## ms      0.28010287 0.002186646 0.2761024 0.003805466 -0.04787150
## mb      0.30255032 0.640630326 0.2837232 0.840650474  0.72126215
##             genizi         car
## md      0.28287857 0.289053925
## richter 0.02561897 0.007137504
## mw      0.11978801 0.086003973
## ms      0.27456078 0.248439249
## mb      0.29715367 0.369365349
## 
## Average coefficients for different model sizes: 
## 
##               1X        2Xs        3Xs        4Xs        5Xs
## md      3.642972  2.5063819  1.7278661  1.2609083  1.0972821
## richter 3.375470 -0.1077359 -0.4982662 -0.5356952 -0.5222478
## mw      5.581534  1.8205076  0.6881174  0.6836889  0.7047378
## ms      3.632778  2.4094065  1.3629567  0.5045514 -0.1807449
## mb      3.702708  3.2636342  2.9579557  2.7766211  2.7010843

plot(importancia)

Ejecutamos diferentes funciones para saber la importancia de los atributos en la tarea de predecir xm.

set.seed(10)

atributos <- c("depth", "year", "month", "lat", "long", "dist")
formula <- as.simple.formula(atributos, "xm")

importancia_SU <- symmetrical.uncertainty(formula, dfMedidas) #medidas para ranking
importancia_IG <- information.gain(formula, dfMedidas)
importancia_GR <- gain.ratio(formula, dfMedidas)
resultado_ranking <- data.frame(importancia_SU, importancia_IG, importancia_GR)
names(resultado_ranking) <- c("Symmetrical Uncertainty", "Information Gain", "Gain Ratio")
round(resultado_ranking,2)

plot(resultado_ranking$`Symmetrical Uncertainty`, resultado_ranking$`Information Gain`)

text(x = resultado_ranking$`Symmetrical Uncertainty`, y = resultado_ranking$`Information Gain`+ 0.0075, cex = 0.75, labels = row.names(resultado_ranking), col="red")

plot(resultado_ranking$`Symmetrical Uncertainty`, resultado_ranking$`Gain Ratio`)

text(x = resultado_ranking$`Symmetrical Uncertainty`, y = resultado_ranking$`Gain Ratio`+
0.0075, cex = 0.75, labels = row.names(resultado_ranking), col="red")

plot(resultado_ranking$`Information Gain`, resultado_ranking$`Gain Ratio`)

text(x = resultado_ranking$`Information Gain`, y = resultado_ranking$`Gain Ratio`+
0.0075, cex = 0.75, labels = row.names(resultado_ranking), col="red")

Además, también usamos random forest.

rf_imp <- random.forest.importance(formula, dfMedidas, importance.type = 1)
print(rf_imp)

##       attr_importance
## depth       118.01025
## year        202.18128
## month        10.57445
## lat          34.17147
## long         79.21607
## dist         17.89221

En todas se puede observar que la variable más significativa es el año, seguido de la profundidad, y después las demás, siendo el mes la menos significativa.

Intentamos encontrar la mejor combinación de atributos, indicándonos CSF que la mejor combinación es la profundidad sin ningún atributo adicional, e indicándonos CON que la mejor combinación es la compuesta por: “depth”, “year”, “lat”, “long” y “dist”. Por tanto, consideramos que no hay un consenso, y que no se puede sacar ninguna conclusión.

subset_CFS <- cfs(formula, dfMedidas)
subset_CON <- consistency(formula, dfMedidas)
print("CFS")

## [1] "CFS"

subset_CFS

## [1] "depth"

print("CON")

## [1] "CON"

subset_CON

## [1] "depth" "year"  "lat"   "long"  "dist"

Llegados a este punto, vamos a aplicar reducción de la dimensionalidad. Para ello vamos a utilizar además de la latitud, longitud y profundidad, “country”, “city”, “area” y “direction”, para ver que tal predicen el intervalo xm de los terremotos. Para poder utilizar estos últimos, necesitamos convertirlos en integer, quitar los valores nulos, y estandarizar todas las variables para que sean comparables.

atr <- c("depth", "country", "city", "area", "direction", "lat", "long")
df <- dfMedidas[ , names(dfMedidas) %in% c(atr, "xm_Interval")]

df$country <- as.integer(df$country)
df$city <- as.integer(df$city)
df$area <- as.integer(df$area)
df$direction <- as.integer(df$direction)
df$xm_Interval <- as.integer(df$xm_Interval)

df <- na.omit(df) #eliminamos instancias si no tienen algún valor

df_scale <- as.data.frame(scale(df)) # estandarizamos los atributos

Ahora que tenemos preparados los datos, vamos a sacar los componentes principales y su importancia.

df.pca <- princomp(df_scale[,-8], cor = TRUE, scores = TRUE, covmat = NULL)

summary(df.pca)

## Importance of components:
##                           Comp.1    Comp.2    Comp.3    Comp.4     Comp.5
## Standard deviation     1.7248537 1.0648893 0.9401792 0.8558231 0.69005497
## Proportion of Variance 0.4250172 0.1619985 0.1262767 0.1046333 0.06802512
## Cumulative Proportion  0.4250172 0.5870157 0.7132924 0.8179257 0.88595080
##                           Comp.6    Comp.7
## Standard deviation     0.6660965 0.5955333
## Proportion of Variance 0.0633835 0.0506657
## Cumulative Proportion  0.9493343 1.0000000

fviz_screeplot(df.pca, addlabels = TRUE)

Se puede ver que el primer componente principal explica el 42.5% de la variabilidad de los datos, y junto al segundo (16.2%), suma el 58.7%.

Vemos el número óptimo de componentes principales a tener en cuenta.

ev <- get_eig(df.pca)[1] #los autovalores
ap <- parallel(subject=nrow(df_scale),var=ncol(df_scale[, -8]), rep=100, cent=.05)
nS <- nScree(x=ev$eigenvalue, aparallel=ap$eigen$qevpea)
plotnScree(nS)

Tres de las cuatro medidas nos indican que lo mejor es tener en cuenta los dos primeros, por lo que las dimensiones se reducen de 7 a 2, realizando combinaciones lineales de los anteriores atributos.

Obtenemos un gráfico biplot, donde vamos a examinar los atributos con su signo, distancia al origen y ángulo con los ejes de los dos primeros componentes principales. De esta forma podremos obtener la correlación y calidad de representación que existe entre los atributos y dichos componentes principales.

Los atributos que estén correlados entre sí, se encontrarán próximos en ángulo, y su calidad de representación, correlación y contribución con cada dimensión, se medirán teniendo en cuenta la proximidad de los atributos a cada eje de coordenadas, además de lo cerca que estén dichos atributos del círculo unidad, siendo el eje X la primera componente principal y el eje Y la segunda.

fviz_pca_var(df.pca)

Se puede apreciar que todos excepto la latitud y la profundidad contribuyen bastante al primer componente principal, aunque la longitud al estar lejos del circulo unidad, no tanto.

En lo referente al segundo componente principal, se puede apreciar que la latitud y la profundidad contribuyen bastante, teneiendo la profundidad una correlación inversa con la latitud, y se podría decir que “country” también contribuye, al estar tan cerca del ciculo uunidad, aunque mucho menos que las otras dos.

Vemos la contribución de los atributos a los dos primeros componentes principales por medio de gráficos.

fviz_contrib(df.pca, choice = "var", axes = 1)

fviz_contrib(df.pca, choice = "var", axes = 2)

Se ratifica la explicación que hemos realizado del biplot.

Ejecutamos otro biplot, pero esta vez con las instancias y los atributos, con diferente color cada intervalo de xm.

fviz_pca_biplot(df.pca, habillage = dfMedidas$xm_Interval)

Se aprecia que las tres clases están completamente solapadas.

Al representarlo con elipses, se sigue apreciando que las tres clases están completamente solapadas.

fviz_pca_ind(df.pca, label="none", habillage = dfMedidas$xm_Interval, addEllipses=TRUE, ellipse.level=0.90)

En vista de que nuestros datos no son muy buenos para hacer clasificación, vamos a probar con la clusterización.

Usamos los mismos atributos que en la reducción de atributos, solo que esta vez utilizamos “xm” en vez de “xm_Interval”, y al igual que antes, pasamos los atributos necesarios a integer, eliminamos los valores nulos, y estandarizamos todos los atributos.

df_CLUS <- dfMedidas[ , names(dfMedidas) %in% c(atr, "xm")]

df_CLUS$country <- as.integer(df_CLUS$country)
df_CLUS$city <- as.integer(df_CLUS$city)
df_CLUS$area <- as.integer(df_CLUS$area)
df_CLUS$direction <- as.integer(df_CLUS$direction)
df_CLUS$xm <- as.integer(df_CLUS$xm)

df_CLUS <- na.omit(df_CLUS) #eliminamos instancias si no tienen algún valor

df_CLUS_scale <- as.data.frame(scale(df_CLUS)) # estandarizamos los atributos

Realizamos una agrupación jerárquica de las instancias de datos.

distancia <- dist(df_CLUS_scale, method = "euclidean") #matriz de distancias euclídeas

clus_hc <- hclust(distancia, method= "ward.D2")

plot(clus_hc)

Se aprecian una gran cantidad de agrupaciones. Lo que habría que hacer ahora es determinar cual es el número óptimo de agrupamientos.

Antes de eso, realizaremos varias representaciones. Una de ellas, será otra agrupación jerárquica, pero esta vez sobre el conjunto de datos transpuesto para ver que atributos están más próximos entre sí.

distanciaT <- dist(t(df_CLUS_scale), method = "euclidean") #matriz de distancias euclídeas

clus_hcT <- hclust(distanciaT, method="ward.D")

fviz_dend(clus_hcT)

Se aprecia una agrupación entre la profundidad y la medida xm, y otra entre las medidas de posición, estándo más próximas entre sí por un lado la latitud y la longitud, y por otro las demás medidas, estándo dentro de estas más próximos los atributos “country” y “direction”, y la agrupación de estos con “city”, para finalmente estar esta agrupación próxima con “area”.

Hacemos una representación de las instancias en los ejes de los dos primeros componentes principales.

fviz_pca_ind(prcomp(df_CLUS_scale), title = "Componentes Principales", geom = "point")

Se puede ver que las instancias no se distribuyen de forma aleatoria, pero no se aprecia ninguna agrupación que sea demasiado clara.

En cuanto a la determinación del número óptimo de agrupamientos, debemos destacar que no nos ha sido posible usar ningún método diferente de “CLARA”, el cual sirve para muestras grandes de datos, ya que usando los demás, como “k-means” o “PAM”, nuestro ordenador no era capaz de ejecutar los métodos o representaciones, bloqueándonos el ordenador en muchas ocasiones, y sacándonos un error indicativo de falta de memeoria RAM en otras. También destacar que con el método “CLARA”, a pesar de darnos un resultado satisfactorio, nos tarda varias horas en ejecutar.

Usándo el método “CLARA”, determinamos un posible número óptimo de agrupamientos, el cual es 3.

fviz_nbclust(df_CLUS_scale, clara, method = "gap_stat")

Volvemos a mostrar la agrupación jerárquica, pero esta vez dividiendo los datos en 3 grupos por medio de recuadros de colores.

plot(clus_hc, cex = 0.6) # plot tree
rect.hclust(clus_hc, k = 3, border = 2:5) # add rectangle

Pensamos que tiene bastante sentido dividir los datos en 3 grupos, ya que en la agrupación jerárquica se aprecian los 3 grupos claramente, aunque también podrían haber sido 2 debido a los dos grandes grupos que se observan.

Representamos en los ejes de los dos primeros componentes principales, la asociación de los atributos a las 3 clases obtenidas, usando k-medias.

km.res <- kmeans(df_CLUS_scale, 3, nstart = 25)
fviz_cluster(km.res, data = df_CLUS_scale, ellipse.type = "convex")

Se aprecia bastante solapamiento de las clases 1 y 2, y un poco de estas con la 3.

Probamos a repetir la representación, pero dividiendo los atributos en 2 clases.

km.res <- kmeans(df_CLUS_scale, 2, nstart = 25)
fviz_cluster(km.res, data = df_CLUS_scale, ellipse.type = "convex")

No se aprecia bien cuanto solapamiento hay, pero se ve que están al menos algo solapadas.

Ahora procedemos a hacer un estudio completo con el método “CLARA”.

Primero veremos que tal se ajustan los datos a las 3 clases que hemos obtenido.

df_clara <- clara(df_CLUS_scale, 3)

Representamos la asociación de los atributos a las 3 clases obtenidas, en los ejes de los dos primeros componentes principales.

fviz_cluster(df_clara, stand = T, geom = "point", pointsize = 1)

Se aprecia solapamiento entre las 3 clases, sobre todo entre la 2 y la 3, pero se observa que hay muchos datos que está claro a que clase pertenecen.

Representamos dos gráficos “Silhouete” diferentes. Si las medidas se acercan a 1, significará que los datos están bien agrupados, y si se acercan a -1 que están mal agrupados.

plot(silhouette(df_clara),  col = 2:3, main = "Silhouette plot")

fviz_silhouette(df_clara)

##   cluster size ave.sil.width
## 1       1   23          0.33
## 2       2   14          0.28
## 3       3    9          0.67

Determinamos que los resultados son bastante buenos, y que los datos están bastante bien agrupados, ya que casi todos los valores de silueta son positivos, incluso llegando a 0.67 el de la clase 3.

Si agruparamos los datos en 2 clases estos serían los resultados:

df_clara2 <- clara(df_CLUS_scale, 2)

Representamos la asociación a 2 clases de los atributos, en los ejes de los dos primeros componentes principales.

fviz_cluster(df_clara2, stand = T, geom = "point", pointsize = 1)

Se sique apreciando solapamiento entre las clases, aunque también se observa en este caso que hay muchos datos que está claro a que clase pertenecen.

Representamos los dos gráficos “Silhouete”.

plot(silhouette(df_clara2),  col = 2:3, main = "Silhouette plot")

fviz_silhouette(df_clara2)

##   cluster size ave.sil.width
## 1       1   19          0.36
## 2       2   25          0.15

Se puede apreciar que las medias son peores que cuando dividíamos los datos en 3 clases, ya que se observan más medidas negativas y la media es inferior siendo esta 0.24 frente a los 0.36 de la situación con 3 clases.

Por tanto consideramos que lo más óptimo es hacer 3 agrupaciones.

Por último, vamos a realizar un estudio acerca de como han cambiado las medidas xm de los terremotos a lo largo de los años de los que tenemos constancia, usando para ello series temporales.

Para llevar acabo esta tarea, hemos tenido que crear un nuevo conjunto de datos, en el que la primera columna contiene los años, y las 12 restantes las medias de xm para cada mes de cada año. (Hemos añadido los meses que no aparecían en nuestro datagrama, y hemos puesto ceros como valor de su media de xm para ese mes y ese año)

dfMediaXm <- read.csv("earthquake_xm_average.csv", sep = ",", dec = ".", header = TRUE)

dfMediaXm

Pensamos que el hecho de que haya meses sin medidas de terremoto se podría deber a dos razones: que no hubo ningún terremoto en esos meses, o que no se tiene constancia de ello. Teniendo en cuenta que esto solo ocurre hasta 1963 y luego en 2017, está bastante claro que se debe a lo segundo.

Antes de empezar con el estudio, decir que no hemos tenido en cuenta los años que tienen meses con valor 0, ya que pensamos que no representan adecuadamente la realidad. Por tanto, nuestro modelo empezará en el año 1964 y terminará en el año 2016, siendo este último el año que intentaremos predecir a partir de los anteriores.

Tomamos el año 2016 como el año a predecir, y creamos un vector con las medidas mensuales de los años restantes.

df_xm_2016 <- dfMediaXm[107,2:13] #Usamos el año 2016 para comparar
df_xm_MED <- as.vector(t(dfMediaXm[55:106,2:13]))

Transformamos los datos a la estructura de serie temporal con ts, y representamos dos gráficos; uno con todas las medidas desde 1964 hasta 2015; y otro con las medidas de 2016.

ts_xm_MED <- ts(df_xm_MED, frequency = 12, start = 1964)
ts_xm_2016 <- ts(t(df_xm_2016), frequency = 12, start = 2016)
ts.plot(ts_xm_MED)

ts.plot(ts_xm_2016)

De momento no podemos sacar ninguna conclusión, pero se aprecia que por lo general a medida que pasan los años la media de xm va disminuyendo de forma muy irregular, hasta que en algún momento entre 1995 y el 2000, empieza a aumentar poco a poco.

Sacamos otros dos gráficos, pero esta vez representando los valores por frecuencia, de manera que cada año sea una serie, con un color más oscuro para los años más antiguos y un color más claro para los más modernos. Uno de ellos está representado en coordenadas lineales, y el otro con una coordenada polar.

ggseasonplot(ts_xm_MED, continuous = TRUE)

ggseasonplot(ts_xm_MED, continuous = TRUE, polar = TRUE)

Se puede apreciar que en los años más antiguos hay bastante más irregularidad en las medidas que en los más modernos, además de apreciarse de forma clara que las medidas de estos últimos años son más bajas que las de los primeros.

Mostramos los valores por subseries mensuales a través de otro gráfico.

ggsubseriesplot(ts_xm_MED)

De este gráfico, se puede resaltar que la media de cada mes es muy parecida, por lo que se puede afirmar que no hay ninguna relación entre los meses y la magnitud de los terremotos. Además, se aprecia claramente lo que expusimos en la primera gráfica.

En el siguiente gráfico mostramos una representación del gráfico de dispersión del dato en un instante t frente al dato retrasado en un instante t-lag(n). Los colores indican los meses.

Con este gráfico intentaremos encontrar cual es la frecuencia con que nuestros datos se repiten o se parecen, si es que hay alguna.

gglagplot(ts_xm_MED, do.lines = FALSE, lags = 12) + geom_point(size = 0.1, stroke = 0)

La frecuencia se puede estimar a partir del valor de retraso (lag) en el que la correlación es máxima, pero como se puede observar, todos tienen una correlación similar, por lo que no se puede afirmar que haya ningún tipo de frecuencia.

Gracias a la función acf, calculamos la autocovarianza y la autocorrelación parcial.

ggAcf(ts_xm_MED)

Sigue sin apreciarse ninguna frecuencia, ya que nuestros datos rara vez se repiten o se parecen.

Representamos la serie en el dominio de las frecuencias para intentar determinar la frecuencia dominante.

spectrum(ts_xm_MED)

El mayor pico ocurre fuera de nuestra frecuencia de 12 meses, por lo que al menos se puede afirmar que los datos no se van a parecer si miramos el mismo mes en diferentes años, lo que ratifica nuestra afirmación de que no existe ninguna relación entre los meses y la magnitud de los terremotos.

Volvemos a usar la función acf, pero utilizando los mayores valores de lag que nos permiten nuestros datos.

ggAcf(ts_xm_MED, lag = 700)

Se aprecia claramente que nuestros datos no se repiten, por lo que se puede afirmar que no existe una frecuencia a la que se repitan o se parezcan nuestros datos, a no ser que esa frecuencia sea 1, o en un caso muy improbable, que sea mucho mayor que el periodo de tiempo del que tenemos datos.

Mostramos el resultado de la descomposición de las componentes de la serie temporal en un modelo aditivo, usando medias móviles y regresión local (LOES), y a partir de ello, establecemos una regresión para la tendencia.

ts_xm_MED_desc <- decompose(ts_xm_MED, type = "additive")
autoplot(ts_xm_MED_desc)

ts_xm_MED_stl <- stl(ts_xm_MED, s.window = "periodic")
autoplot(ts_xm_MED_stl)

autoplot(ts_xm_MED_stl$time.series[,2]) +
stat_smooth(method = "lm", aes(color = "Lineal"), se = TRUE, level = 0.95) +
stat_smooth(method = "loess" , aes(color = "LOESS"), se = TRUE, level = 0.95)

Con la regresión lineal se aprecia claramente la disminución de xm a medida que pasan los años, y con la regresión local se observa la tendencia que conseguimos apreciar en un principio, en la que xm disminuye hasta mitad de los años 90, y entonces empieza a aumentar poco a poco.

Ahora vamos a apilcar unos cuantos métodos para intentar predecir las medidas xm del año 2016.

Primero aplicamos un método naive, que predice la parte de tendencia para posteriormente aplicarla a los datos desestacionalizados. Sólo mostramos desde el año 2010, para que se vean bien los datos predichos.

fit_naive_ts <- stlf(ts_xm_MED, method = "naive", h = 12, s.window = 2, robust = TRUE)
autoplot(fit_naive_ts, xlim = c(2010, 2017))

## Scale for 'x' is already present. Adding another scale for 'x', which
## will replace the existing scale.

## Warning: Removed 552 rows containing missing values (geom_path).

Usamos el método Holt-Winters para hacer suavizado exponencial, que da menos peso a los datos cuanto más antiguos sean.

fit_HW_ts <- hw(ts_xm_MED, h = 12)
autoplot(fit_HW_ts, xlim = c(2010, 2017))

## Scale for 'x' is already present. Adding another scale for 'x', which
## will replace the existing scale.

## Warning: Removed 552 rows containing missing values (geom_path).

Los métodos de suavizado exponencial demandandan, que los errores en el intervalo de predicción, estén distribuidos según una normal de media cero, por lo que los errores no están correlados y la varianza es constante.

Se puede utilizar la correlación del error para mejorar la estimación del modelo de predicción usando ARIMA, pero para ello, los datos de la serie se deben preparar.

En las gráficas siguientes mostramos la serie preparada para aplicar ARIMA. En la primera se muestra la transformación de la serie para poder aplicar dicho método, y en la segunda la autocorrelación de la serie transformada.

ts_xm_MED_ARIMA <- diff(ts_xm_MED, differences = 12)
autoplot(ts_xm_MED_ARIMA)

ggAcf(ts_xm_MED_ARIMA)

Ya podemos aplicar ARIMA.

fit_ARIMA_model <- auto.arima(ts_xm_MED, stepwise = FALSE, approximation = FALSE)
fit_ARIMA_ts <- forecast(fit_ARIMA_model)
autoplot(fit_ARIMA_ts, xlim = c(2010, 2017))

## Scale for 'x' is already present. Adding another scale for 'x', which
## will replace the existing scale.

## Warning: Removed 552 rows containing missing values (geom_path).

Finalmente, aplicamos el modelo TBATS, que es una extensión del modelo de suavizado exponencial, y también obtiene un modelo ARIMA.

fit_TBATS_model <- auto.arima(ts_xm_MED)
fit_TBATS_ts <- forecast(fit_TBATS_model)
autoplot(fit_TBATS_ts, xlim = c(2010, 2017))

## Scale for 'x' is already present. Adding another scale for 'x', which
## will replace the existing scale.

## Warning: Removed 552 rows containing missing values (geom_path).

representamos en una gráfica con los valores reales y los predichos por cada uno de los métodos a ver cual predice mejor la realidad.

ts.plot( fit_naive_ts$mean, fit_naive_ts$lower, fit_naive_ts$upper,
fit_HW_ts$mean, fit_HW_ts$lower, fit_HW_ts$upper,
fit_ARIMA_ts$mean, fit_ARIMA_ts$lower, fit_ARIMA_ts$upper,
fit_TBATS_ts$mean, fit_TBATS_ts$lower, fit_TBATS_ts$upper,
ts_xm_2016,
col = c(2,2,0,2,0,3,3, 0, 3, 0, 4, 4, 0, 4, 0, 5, 5, 0, 5, 0, 1),
lty = c(1,2,2,2,2, 1, 2, 2, 2, 2, 1, 2, 2, 2, 2, 1, 2, 2, 2, 2, 1),
lwd = c(2,0.5,0.5,0.5, 0.5, 2, 0.5, 0.5, 0.5, 0.5, 2, 0.5, 0.5, 0.5, 0.5, 2, 0.5, 0.5, 0.5, 0.5, 2),
gpars=list(xaxt="n"),
xlab = "Mes", ylab = "xm medio")
points(seq(2015,2016-(1/12),by=1/12),fit_naive_ts$mean[1:12],pch=19, col = 2)
points(seq(2016,2017-(1/12),by=1/12),fit_HW_ts$mean[1:12],pch= 19, col = 3)
points(seq(2016,2017-(1/12),by=1/12),fit_ARIMA_ts$mean[1:12],pch= 19, col = 4)
points(seq(2016,2017-(1/12), by=1/12),fit_TBATS_ts$mean[1:12],pch= 19,col = 5)
points(seq(2016,2017-(1/12), by=1/12), ts_xm_2016[1:12], pch = 19, col = 1)
axis(1 ,at=seq(2016,2017-(1/12), by=1/12), labels=c("ENE", "FEB", "MAR", "ABR", "MAY",
"JUN", "JUL", "AGO", "SEP", "OCT", "NOV", "DIC"))
legend("topleft", legend = c("Naive","Holt-Winters","ARIMA", "TBATS", "Real"),
col = c(2, 3, 4, 5, 1), lty = 1, text.width = 0.2, cex = 0.5)
title("Predicción media xm 2016", cex.main = 1, sub = "Intervalo 95%", cex.sub = 1)

Se puede ver que ningún modelo se acerca mucho a la realidad, ya que todos menos Naive predicen que las medidas serán mas o menos constantes, cuando en realidad no lo son, y Naive predice por lo general medidas mas bajas que las reales. Pero hay que tener en cuenta que los datos de las medias de xm están guardados con muchos decimales, por lo que la estimación de todos los modelos no se aleja tanto de la realidad.

Lo veremos comparamos la precisión, tanto en entrenamiento como en validación, de las diferentes técnicas utilizadas.

rmse_nv <- forecast::accuracy(fit_naive_ts)[2]
rmse_hw <- forecast::accuracy(fit_HW_ts)[2]
rmse_ar <- forecast::accuracy(fit_ARIMA_model)[2]
rmse_tb <- forecast::accuracy(fit_TBATS_model)[2]
par(mfrow=c(2,1))
barplot(c(rmse_nv, rmse_hw, rmse_ar, rmse_tb), names.arg = c("Naive", "HW", "Arima",
"TBATS"), ylim = c(0, 3), main = "RMSE entrenamiento 1912-2015")
text(c(0.7,1.9,3.1,4.3), 0.3 + round(c(rmse_nv, rmse_hw, rmse_ar, rmse_tb), 1), labels
= round(c(rmse_nv, rmse_hw, rmse_ar, rmse_tb), 1))
barplot(c(rmse(ts_xm_2016, fit_naive_ts$mean), rmse(ts_xm_2016, fit_HW_ts$mean), rmse(ts_xm_2016,
fit_ARIMA_ts$mean), rmse(ts_xm_2016, fit_TBATS_ts$mean)), names.arg = c("Naive", "HW",
"Arima", "TBATS"), ylim = c(0, 3), main = "RMSE año 2016")
text(c(0.7,1.9,3.1,4.3), 0.3 + round(c(rmse(ts_xm_2016, fit_naive_ts$mean), rmse(ts_xm_2016,
fit_HW_ts$mean), rmse(ts_xm_2016, fit_ARIMA_ts$mean), rmse(ts_xm_2016, fit_TBATS_ts$mean)),
1), labels = round(c( rmse(ts_xm_2016, fit_naive_ts$mean), rmse(ts_xm_2016, fit_HW_ts$mean),
rmse(ts_xm_2016, fit_ARIMA_ts$mean), rmse(ts_xm_2016, fit_TBATS_ts$mean)), 1))

Se puede observar que todos excepto Naive en validación (0.2), tienen un error de 0.1, por lo que consideramos que se ha hecho una buena estimación de la media de xm para los meses del año 2016.

Por último, observamos si hay algún dato anómalo.

outliers <- bfast(ts_xm_MED, h=0.15, max.iter=1)
outliers

## 
##   TREND BREAKPOINTS
##   Confidence intervals for breakpoints
##   of optimal 3-segment partition: 
## 
## Call:
## confint.breakpointsfull(object = bp.Vt, het.err = FALSE)
## 
## Breakpoints at observation number:
##   2.5 % breakpoints 97.5 %
## 1   132         133    134
## 2   345         365    367
## 
## Corresponding to breakdates:
##   2.5 %    breakpoints 97.5 % 
## 1 1974(12) 1975(1)     1975(2)
## 2 1992(9)  1994(5)     1994(7)
## 
##   SEASONAL BREAKPOINTS:  None

plot(outliers)

No lo hay, pero se aprecian dos discontinuidades en la tendencia: una alrededor de la mitad de los años 70, donde empieza a disminuir xm más despacio que hasta ese momento, u otra a mitad de los años 90, donde xm empieza a aumentar.

Cómo conclusión del estudio total de la práctica, pensamos que es complicado predecir cuando y donde se producirán los terremotos y de que escala serán, ya que ningún atributo por si sólo ha sido capaz de predecir muy bien estos datos, y hemos descubierto además que la mejor manera de predecir su magnitud es estudiar la magnitud de los terremotos que han acontecido en el pasado.